########################################################################################################################
# PROJECT:    Determinants of COVID-19 Vaccine Fatigue
# OBJECTIVE:  Analysis
# AUTHOR:     TS, JP, VR
# DATE:       2022-08-18 to 2023-01-30
########################################################################################################################

if(!suppressWarnings(require(conjoint))) install.packages("conjoint") ; library(conjoint)
if(!suppressWarnings(require(haven))) install.packages("haven") ; library(haven)
if(!suppressWarnings(require(readxl))) install.packages("readxl") ; library(readxl)
if(!suppressWarnings(require(dplyr))) install.packages("dplyr") ; library(dplyr)
if(!suppressWarnings(require(reshape2))) install.packages("reshape2") ; library(reshape2)
if(!suppressWarnings(require(cregg))) install.packages("cregg") ; library(cregg)
if(!suppressWarnings(require(haven))) install.packages("haven") ; library(haven)
if(!suppressWarnings(require(ggplot2))) install.packages("ggplot2") ; library(ggplot2)
if(!suppressWarnings(require(ggpubr))) install.packages("ggpubr") ; library(ggpubr)
if(!suppressWarnings(require(tidyverse))) install.packages("tidyverse") ; library(tidyverse)

############################## Conjoint analysis ##############################
### Choices - 1. Experiment
conjoint2 <- read_excel("E:/COVID/Conjoint2/Data_AT_Pub.xlsx")
conjoint3 <- read_excel("E:/COVID/Conjoint2/Data_IT_Pub.xlsx")
conjoint2 <- rbind(conjoint2, conjoint3)
conjoint2$Vaccinated[conjoint2$Vaccinated == "Three_or_more"] <- "3+"
conjoint2$Vaccinated1[conjoint2$Vaccinated1 == "Three_or_more"] <- "3+"
conjoint2$Vaccinated[conjoint2$Vaccinated == "One_or_two"] <- "1_or_2"
conjoint2$Vaccinated1[conjoint2$Vaccinated1 == "One_or_two"] <- "1_or_2"
conjoint2$Vaccinated[conjoint2$Vaccinated == "Not"] <- "0"
conjoint2$Vaccinated1[conjoint2$Vaccinated1 == "Not"] <- "0"
data1 <- select(conjoint2, U0, Ex1A, country, Vaccinated, Gender1, Q3, Attitude, Q1)
data2 <- select(conjoint2, U0, Ex1B, country, Vaccinated, Gender2, Q3, Attitude, Q1)
data3 <- select(conjoint2, U0, Ex1A0, country, Vaccinated, Gender3, Q3, Attitude, Q1)
data4 <- select(conjoint2, U0, Ex1B0, country, Vaccinated, Gender4, Q3, Attitude, Q1)
data1 <- dplyr::rename(data1, ID = U0, Profile = Ex1A, Gender = Gender1, Age_group = Q3, GenderP = Q1, Country=country)
data2 <- dplyr::rename(data2, ID = U0, Profile = Ex1B, Gender = Gender2, Age_group = Q3, GenderP = Q1, Country=country)
data3 <- dplyr::rename(data3, ID = U0, Profile = Ex1A0, Gender = Gender3, Age_group = Q3, GenderP = Q1, Country=country)
data4 <- dplyr::rename(data4, ID = U0, Profile = Ex1B0, Gender = Gender4, Age_group = Q3, GenderP = Q1, Country=country)
data5 <- rbind(data1, data2)
data5$Rating <- 1
data5a <- rbind(data3, data4)
data5a$Rating <- 0
data6 <- rbind(data5, data5a)
profile <- read_excel("E:/COVID/Conjoint2/Experiment1.xlsx")
data9 <- merge(data6, profile, by ="Profile")
data9$Variant <- as.factor(data9$Variant)
data9$Vacc <- as.factor(data9$Vacc)
data9$Omic <- as.factor(data9$Omic)
data9$Incen <- as.factor(data9$Incen)
data9$Target <- as.factor(data9$Target)
data9$Motiv <- as.factor(data9$Motiv)
data9$Country <- as.factor(data9$Country)
data9$Vaccinated <- as.factor(data9$Vaccinated)
data9$Gender <- as.factor(data9$Gender)
data9$Age_group <- as.factor(data9$Age_group)
data9$Attitude <- as.factor(data9$Attitude)
data9$GenderP <- as.factor(data9$GenderP)
data9 <- dplyr::rename(data9, Virus_variant = Variant, Vaccines = Vacc, Omicron_adapted = Omic, Incentives = Incen, Motivation=Motiv)
data9$Omicron_adapted <- factor(data9$Omicron_adapted, levels = c("not_adap","adap"))
levels(data9$Omicron_adapted) <- list(Not_adapted  = "not_adap", Adapted = "adap")
data9$Virus_variant <- factor(data9$Virus_variant, levels = c("No change","Escalation", "Decline"))
levels(data9$Virus_variant) <- list(No_change  = "No change", Escalation = "Escalation", Decline = "Decline")
data9$Incentives <- factor(data9$Incentives, levels = c("Free", "Not_free", "Cash", "Voucher"))
levels(data9$Incentives) <- list(Free = "Free", "20_Euro_fee" = "Not_free", "500_Euro_cash"  = "Cash", "500_Euro_voucher"  = "Voucher")

data9$Motiv <- factor(data9$Motivation, levels = c("Infect_Ego",  "Econ_Ego","Serious_Ego", "Community", "Serious_Friend", "Efficacy_Ego", "Econ_Socio", "Health_System"))
levels(data9$Motivation) <- list("Risk_re-infection" = "Infect_Ego", "Risk_unable_to_work" = "Econ_Ego", "Risk_severe_disease" = "Serious_Ego", "Community_spirit" = "Community", "Protect_friends" = "Serious_Friend", "Self_efficacy" = "Efficacy_Ego", Risk_lockdown = "Econ_Socio", Protect_health_system = "Health_System")
f1 <- Rating ~ Virus_variant + Vaccines + Omicron_adapted + Incentives + Motivation
mm <- cj(na.omit(data9), f1, id = ~ID, estimate = "mm", by = ~Country)
a <- plot(mm, group="Country", vline = 0.5)
a
data9$Country <- factor(data9$Country, levels = c("AT","IT"))
amces <- cj(data9, f1, id = ~ID, by = ~Country)
p <- plot(amces, group = "Country", xlim = c(-0.25, 0.25))
p1 <- p + scale_color_manual(values=c("#d7191c", "#2c7bb6"), na.translate = F)
p1
write.table(amces, "E:/COVID/Conjoint2/amce1_Choice_both.txt", sep="\t")

#AMCEs by Vaccination Status
amces <- cj(data9, f1, id = ~ID, by = ~Vaccinated)
p2 <- plot(amces, group="Vaccinated", xlim = c(-0.25, 0.25))
p2
write.table(amces, "E:/COVID/Conjoint2/amce1_Choice_both.txt", sep="\t")

#AMCEs by Gender
data9G <- subset(data9, (GenderP == "female") | (GenderP == "male"))
data9G$GenderP <- factor(data9G$GenderP, levels = c("female","male"))
amces <- cj(data9G, f1, id = ~ID, by = ~GenderP)
pG1 <- plot(amces, group="GenderP", xlim = c(-0.25, 0.25))
pG1 <- pG1 + scale_color_manual(values=c("pink", "darkblue"), na.translate = F) 
pG1 <- pG1 + ggtitle("Evaluation of campaign (binary choice)") +
  theme(plot.title = element_text(hjust = 0.5))
pG1
write.table(amces, "E:/COVID/Conjoint2/amce1_Choice_both.txt", sep="\t")

#Has Gender an effect?
f1 <- Rating ~ Virus_variant + Vaccines + Omicron_adapted + Incentives + Motivation + Gender
mm <- cj(na.omit(data9), f1, id = ~ID, estimate = "mm", by = ~Country)
a <- plot(mm, group="Country", vline = 0.5)
a
amces <- cj(data9, f1, id = ~ID, by = ~Country)
p <- plot(amces, group="Country", xlim = c(-0.5, 0.5))
p

### Choices - 2. Experiment
data1 <- select(conjoint2, U0, Ex2A, country, Vaccinated, Gender1, Q1)
data2 <- select(conjoint2, U0, Ex2B, country, Vaccinated, Gender2, Q1)
data3 <- select(conjoint2, U0, Ex2A0, country, Vaccinated, Gender3, Q1)
data4 <- select(conjoint2, U0, Ex2B0, country, Vaccinated, Gender4, Q1)
data1 <- dplyr::rename(data1, ID = U0, Profile = Ex2A, Country=country, Gender = Gender1, GenderP = Q1)
data2 <- dplyr::rename(data2, ID = U0, Profile = Ex2B, Country=country, Gender = Gender2, GenderP = Q1)
data3 <- dplyr::rename(data3, ID = U0, Profile = Ex2A0, Country=country, Gender = Gender3, GenderP = Q1)
data4 <- dplyr::rename(data4, ID = U0, Profile = Ex2B0, Country=country, Gender = Gender4, GenderP = Q1)
data5 <- rbind(data1, data2)
data5$Rating <- 1
data5a <- rbind(data3, data4)
data5a$Rating <- 0
data6 <- rbind(data5, data5a)
profile <- read_excel("E:/COVID/Conjoint2/Experiment2.xlsx")
data9 <- merge(data6, profile, by ="Profile")
data9$Cons <- as.factor(data9$Cons)
data9$Celeb <- as.factor(data9$Celeb)
data9$LongCo <- as.factor(data9$LongCo)
data9$GreenPa <- as.factor(data9$GreenPa)
data9$Mand <- as.factor(data9$Mand)
data9$Country <- as.factor(data9$Country)
data9$Vaccinated <- as.factor(data9$Vaccinated)
data9$Gender <- as.factor(data9$Gender)
data9$GenderP <- as.factor(data9$GenderP)
data9$Cons <- factor(data9$Cons, levels = c("Cons_scie","Cons_phys","False_bal_scie", "False_bal_phys"))
data9$Celeb <- factor(data9$Celeb, levels = c("Celeb_not_vacc","Celeb_vacc","Celeb_waits", "Celeb_ill"))
data9$LongCo <- factor(data9$LongCo, levels = c("1Percent","5Percent","10Percent", "20Percent"))
data9$GreenPa <- factor(data9$GreenPa, levels = c("Not_need","Needed"))
data9$Mand <- factor(data9$Mand, levels = c("None","50 Years-100 Euro", "18 Years-1500 Euro"))

data9 <- dplyr::rename(data9, Consensus = Cons, Celebrity = Celeb, Long_COVID = LongCo, Green_pass = GreenPa, Vaccine_mandate = Mand)
levels(data9$Consensus) <- list("Consensus_scientists" = "Cons_scie", "Consensus_physicians" = "Cons_phys", "Dissensus_scientists" = "False_bal_scie", "Dissensus_physicians" = "False_bal_phys")
levels(data9$Celebrity) <- list("Refuses_vaccine" = "Celeb_not_vacc", "Endorses_vaccine" = "Celeb_vacc", "Waits_for_new_vaccine" = "Celeb_waits", "Infected_regrets_hesitancy" = "Celeb_ill")
levels(data9$Long_COVID) <- list("1_percent" = "1Percent", "5_percent" = "5Percent", "10_percent" ="10Percent", "20_percent" = "20Percent")
levels(data9$Vaccine_mandate) <- list("No" = "None", "Age_50+_fine_100Euro" = "50 Years-100 Euro", "Age_18+_fine_1500Euro" = "18 Years-1500 Euro")
f1 <- Rating ~ Consensus + Celebrity + Long_COVID + Green_pass + Vaccine_mandate

mm <- cj(data9, f1, id = ~ID, estimate = "mm", by = ~Country)
a <- plot(mm, group="Country", vline = 0.5)
a
amces <- cj(data9, f1, id = ~ID, by = ~Country)
p5 <- plot(amces, group="Country", xlim = c(-0.25, 0.25))
p5 <- p5 + scale_color_manual(values=c("#d7191c", "#2c7bb6"), na.translate = F)
p5

# AMCEs by Vaccination Status
amces <- cj(data9, f1, id = ~ID, by = ~Vaccinated)
p6 <- plot(amces, group="Vaccinated", xlim = c(-0.25, 0.25))
p6

# AMCEs by Gender
data9G <- subset(data9, (GenderP == "female") | (GenderP == "male"))
data9G$GenderP <- factor(data9G$GenderP, levels = c("female","male"))
amces <- cj(data9G, f1, id = ~ID, by = ~GenderP)
pG2 <- plot(amces, group="GenderP", xlim = c(-0.25, 0.25))
pG2 <- pG2 + scale_color_manual(values=c("pink", "darkblue"), na.translate = F)
pG2 <- pG2 + ggtitle("Trust in vaccine (binary choice)") +
  theme(plot.title = element_text(hjust = 0.5))
pG2
write.table(amces, "E:/COVID/Conjoint2/amce2_Choice_both.txt", sep="\t")

# Has Gender an effect?
f1 <- Rating ~ Consensus + Celebrity + Long_COVID + Green_pass + Vaccine_mandate + Gender
mm <- cj(na.omit(data9), f1, id = ~ID, estimate = "mm", by = ~Country)
a <- plot(mm, group="Country", vline = 0.5)
a
amces <- cj(data9, f1, id = ~ID, by = ~Country)
p <- plot(amces, group="Country", xlim = c(-0.5, 0.5))
p

### Ratings - 1. Experiment
conjoint2 <- read_excel("E:/COVID/Conjoint2/Data_AT_Pub.xlsx")
conjoint3 <- read_excel("E:/COVID/Conjoint2/Data_IT_Pub.xlsx")
conjoint2 <- rbind(conjoint2, conjoint3)
conjoint2$Vaccinated[conjoint2$Vaccinated == "Three_or_more"] <- "3+"
conjoint2$Vaccinated1[conjoint2$Vaccinated1 == "Three_or_more"] <- "3+"
conjoint2$Vaccinated[conjoint2$Vaccinated == "One_or_two"] <- "1_or_2"
conjoint2$Vaccinated1[conjoint2$Vaccinated1 == "One_or_two"] <- "1_or_2"
conjoint2$Vaccinated[conjoint2$Vaccinated == "Not"] <- "0"
conjoint2$Vaccinated1[conjoint2$Vaccinated1 == "Not"] <- "0"
data1 <- select(conjoint2, U0, Q18A1, Q20A1, country, Vaccinated, Gender1, Q1)
data2 <- select(conjoint2, U0, Q18A2, Q20A2, country, Vaccinated, Gender2, Q1)
data3 <- select(conjoint2, U0, Q21A1, Q23A1, country, Vaccinated, Gender3, Q1)
data4 <- select(conjoint2, U0, Q21A2, Q23A2, country, Vaccinated, Gender4, Q1)

data1 <- dplyr::rename(data1, ID = U0, Profile = Q18A1, Rating = Q20A1, Country=country, Gender = Gender1, GenderP = Q1)
data2 <- dplyr::rename(data2, ID = U0, Profile = Q18A2, Rating = Q20A2, Country=country, Gender = Gender2, GenderP = Q1)
data3 <- dplyr::rename(data3, ID = U0, Profile = Q21A1, Rating = Q23A1, Country=country, Gender = Gender3, GenderP = Q1)
data4 <- dplyr::rename(data4, ID = U0, Profile = Q21A2, Rating = Q23A2, Country=country, Gender = Gender4, GenderP = Q1)
data5 <- rbind(data1, data2, data3, data4)

profile <- read_excel("E:/COVID/Conjoint2/Experiment1.xlsx")
data9 <- merge(data5, profile, by ="Profile")
data9$Variant <- as.factor(data9$Variant)
data9$Vacc <- as.factor(data9$Vacc)
data9$Omic <- as.factor(data9$Omic)
data9$Incen <- as.factor(data9$Incen)
data9$Target <- as.factor(data9$Target)
data9$Motiv <- as.factor(data9$Motiv)
data9$Country <- as.factor(data9$Country)
data9$Vaccinated <- as.factor(data9$Vaccinated)
data9$Gender <- as.factor(data9$Gender)
data9$GenderP <- as.factor(data9$GenderP)
data9$Variant <- factor(data9$Variant, levels = c("No change","Escalation", "Decline"))
data9$Omic <- factor(data9$Omic, levels = c("not_adap","adap"))
data9$Incen <- factor(data9$Incen, levels = c("Free", "Not_free", "Cash", "Voucher"))
data9$Motiv <- factor(data9$Motiv, levels = c("Infect_Ego",  "Econ_Ego","Serious_Ego", "Community", "Serious_Friend", "Efficacy_Ego", "Econ_Socio", "Health_System"))

data9 <- dplyr::rename(data9, Virus_variant = Variant, Vaccines = Vacc, Omicron_adapted = Omic, Incentives = Incen, Motivation=Motiv)
data9$Omicron_adapted <- factor(data9$Omicron_adapted, levels = c("not_adap","adap"))
levels(data9$Omicron_adapted) <- list(Not_adapted  = "not_adap", Adapted = "adap")
data9$Virus_variant <- factor(data9$Virus_variant, levels = c("No change","Escalation", "Decline"))
levels(data9$Virus_variant) <- list(No_change  = "No change", Escalation = "Escalation", Decline = "Decline")
data9$Incentives <- factor(data9$Incentives, levels = c("Free", "Not_free", "Cash", "Voucher"))
levels(data9$Incentives) <- list(Free = "Free", "20_Euro_fee" = "Not_free", "500_Euro_cash"  = "Cash", "500_Euro_voucher"  = "Voucher")

data9$Motiv <- factor(data9$Motivation, levels = c("Infect_Ego",  "Econ_Ego","Serious_Ego", "Community", "Serious_Friend", "Efficacy_Ego", "Econ_Socio", "Health_System"))
levels(data9$Motivation) <- list("Risk_re-infection" = "Infect_Ego", "Risk_unable_to_work" = "Econ_Ego", "Risk_severe_disease" = "Serious_Ego", "Community_spirit" = "Community", "Protect_friends" = "Serious_Friend", "Self_efficacy" = "Efficacy_Ego", Risk_lockdown = "Econ_Socio", Protect_health_system = "Health_System")
f1 <- Rating ~ Virus_variant + Vaccines + Omicron_adapted + Incentives + Motivation
y <- plot(mm(data9, f1, id = ~ID), vline = 0.5)
y
mm <- cj(na.omit(data9), f1, id = ~ID, estimate = "mm", by = ~Country)
a <- plot(mm, group="Country")
a
amces <- cj(na.omit(data9), f1, id = ~ID, by = ~Country)
p3 <- plot(amces, group="Country", xlim = c(-1, 1))
p3 <- p3 + scale_color_manual(values=c("#d7191c", "#2c7bb6"), na.translate = F)
write.table(amces, "E:/COVID/Conjoint2/amce1_both.txt", sep="\t")
p3

# AMCEs by Vaccination Status
amces <- cj(data9, f1, id = ~ID, by = ~Vaccinated)
p4 <- plot(amces, group="Vaccinated", xlim = c(-1, 1))
p4

# AMCEs by Gender
data9G <- subset(data9, (GenderP == "female") | (GenderP == "male"))
data9G$GenderP <- factor(data9G$GenderP, levels = c("female","male"))
amces <- cj(data9G, f1, id = ~ID, by = ~GenderP)
pG3 <- plot(amces, group="GenderP", xlim = c(-1, 1))
pG3 <- pG3 + scale_color_manual(values=c("pink", "darkblue"), na.translate = F)
pG3 <- pG3 + ggtitle("Likelihood to get vaccinated (ratings)") +
  theme(plot.title = element_text(hjust = 0.5))
pG3
write.table(amces, "E:/COVID/Conjoint2/amce2_Choice_both.txt", sep="\t")

figure1 <- ggarrange(p1, p3, p2, p4, labels = c("A", "B", "C", "D"), ncol = 2, nrow = 2)
figure1
annotate_figure(figure1, top = text_grob(" Figure 1.                      Evaluation of campaign (binary choice)                                    Likelihood to get vaccinated (ratings)    ", , color = "black", face = "bold", size = 12), left = text_grob("                 By vaccination status                                                                                                                    By country", color = "black", face = "bold", rot = 90), bottom = text_grob("AMCEs = Average Marginal Component Effects     ", color = "black", hjust = 1, x = 1, face = "italic", size = 10))

### Ratings - 2. Experiment
data1 <- select(conjoint2, U0, Q27A1, Q29A1, country, Vaccinated, Gender1, Q1)
data2 <- select(conjoint2, U0, Q27A2, Q29A2, country, Vaccinated, Gender2, Q1)
data3 <- select(conjoint2, U0, Q30A1, Q32A1, country, Vaccinated, Gender3, Q1)
data4 <- select(conjoint2, U0, Q30A2, Q32A2, country, Vaccinated, Gender4, Q1)
data1 <- dplyr::rename(data1, ID = U0, Profile = Q27A1, Rating = Q29A1, Country=country, Gender = Gender1, GenderP = Q1)
data2 <- dplyr::rename(data2, ID = U0, Profile = Q27A2, Rating = Q29A2, Country=country, Gender = Gender2, GenderP = Q1)
data3 <- dplyr::rename(data3, ID = U0, Profile = Q30A1, Rating = Q32A1, Country=country, Gender = Gender3, GenderP = Q1)
data4 <- dplyr::rename(data4, ID = U0, Profile = Q30A2, Rating = Q32A2, Country=country, Gender = Gender4, GenderP = Q1)
data5 <- rbind(data1, data2, data3, data4)
profile <- read_excel("E:/COVID/Conjoint2/Experiment2.xlsx")
data9 <- merge(data5, profile, by ="Profile")
data9$Cons <- as.factor(data9$Cons)
data9$Celeb <- as.factor(data9$Celeb)
data9$LongCo <- as.factor(data9$LongCo)
data9$GreenPa <- as.factor(data9$GreenPa)
data9$Mand <- as.factor(data9$Mand)
data9$Country <- as.factor(data9$Country)
data9$Vaccinated <- as.factor(data9$Vaccinated)
data9$Gender <- as.factor(data9$Gender)
data9$GenderP <- as.factor(data9$GenderP)
data9$Cons <- factor(data9$Cons, levels = c("Cons_scie","Cons_phys","False_bal_scie", "False_bal_phys"))
data9$Celeb <- factor(data9$Celeb, levels = c("Celeb_not_vacc","Celeb_vacc","Celeb_waits", "Celeb_ill"))
data9$LongCo <- factor(data9$LongCo, levels = c("1Percent","5Percent","10Percent", "20Percent"))
data9$GreenPa <- factor(data9$GreenPa, levels = c("Not_need","Needed"))
data9$Mand <- factor(data9$Mand, levels = c("None","50 Years-100 Euro", "18 Years-1500 Euro"))

data9 <- dplyr::rename(data9, Consensus = Cons, Celebrity = Celeb, Long_COVID = LongCo, Green_pass = GreenPa, Vaccine_mandate = Mand)
levels(data9$Consensus) <- list("Consensus_scientists" = "Cons_scie", "Consensus_physicians" = "Cons_phys", "Dissensus_scientists" = "False_bal_scie", "Dissensus_physicians" = "False_bal_phys")
levels(data9$Celebrity) <- list("Refuses_vaccine" = "Celeb_not_vacc", "Endorses_vaccine" = "Celeb_vacc", "Waits_for_new_vaccine" = "Celeb_waits", "Infected_regrets_hesitancy" = "Celeb_ill")
levels(data9$Long_COVID) <- list("1_percent" = "1Percent", "5_percent" = "5Percent", "10_percent" ="10Percent", "20_percent" = "20Percent")
levels(data9$Vaccine_mandate) <- list("No" = "None", "Age_50+_fine_100Euro" = "50 Years-100 Euro", "Age_18+_fine_1500Euro" = "18 Years-1500 Euro")
f1 <- Rating ~ Consensus + Celebrity + Long_COVID + Green_pass + Vaccine_mandate

y <- plot(mm(data9, f1, id = ~ID), vline = 0.5)
y
mm <- cj(na.omit(data9), f1, id = ~ID, estimate = "mm", by = ~Country)
a <- plot(mm, group="Country")
a
amces <- cj(na.omit(data9), f1, id = ~ID, by = ~Country)
p7 <- plot(amces, group="Country", xlim=c(-0.6, 0.6))
p7 <- p7 + scale_color_manual(values=c("#d7191c", "#2c7bb6"), na.translate = F)
p7
write.table(amces, "E:/COVID/Conjoint2/amce1_both.txt", sep="\t")

# AMCEs by Vaccination Status
amces <- cj(data9, f1, id = ~ID, by = ~Vaccinated)
p8 <- plot(amces, group="Vaccinated", xlim = c(-0.6, 0.6))
p8

# AMCEs by Gender
data9G <- subset(data9, (GenderP == "female") | (GenderP == "male"))
data9G$GenderP <- factor(data9G$GenderP, levels = c("female","male"))
amces <- cj(data9G, f1, id = ~ID, by = ~GenderP)
pG4 <- plot(amces, group="GenderP", xlim = c(-1, 1))
pG4 <- pG4 + scale_color_manual(values=c("pink", "darkblue"), na.translate = F)
pG4 <- pG4 + ggtitle("Likelihood to get vaccinated (binary choice)") +
  theme(plot.title = element_text(hjust = 0.5))
pG4
write.table(amces, "E:/COVID/Conjoint2/amce2_Choice_both.txt", sep="\t")

figure2 <- ggarrange(p5, p7, p6, p8, labels = c("A", "B", "C", "D"), ncol = 2, nrow = 2)
figure2
annotate_figure(figure2, top = text_grob("Figure 2.                               Trust in vaccine (binary choice)                                           Likelihood to get vaccinated (ratings)    ", , color = "black", face = "bold", size = 12), left = text_grob("       By vaccination status                                                                                                                     By country", color = "black", face = "bold", rot = 90), bottom = text_grob("AMCEs = Average Marginal Component Effects     ", color = "black", hjust = 1, x = 1, face = "italic", size = 10))

figure1G <- ggarrange(pG1, pG3, labels = c("A", "B"), ncol = 2, nrow = 1)
figure1G
figure2G <- ggarrange(pG2, pG4, labels = c("A", "B"), ncol = 2, nrow = 1)
figure2G

############################## Boxplots ##############################
### Boxplot Trust
data10 <- select (conjoint2, Q15A1, Q15A2, Q15A3, Q15A4, Q15A5, Q15A6, Q15A7, Vaccinated, country)
data10$Q15A1[data10$Q15A1 == "88"] <- NA
data10$Q15A1[data10$Q15A1 == "99"] <- NA
data10$Q15A2[data10$Q15A2 == "88"] <- NA
data10$Q15A2[data10$Q15A2 == "99"] <- NA
data10$Q15A3[data10$Q15A3 == "88"] <- NA
data10$Q15A3[data10$Q15A3 == "99"] <- NA
data10$Q15A4[data10$Q15A4 == "88"] <- NA
data10$Q15A4[data10$Q15A4 == "99"] <- NA
data10$Q15A5[data10$Q15A5 == "88"] <- NA
data10$Q15A5[data10$Q15A5 == "99"] <- NA
data10$Q15A6[data10$Q15A6 == "88"] <- NA
data10$Q15A6[data10$Q15A6 == "99"] <- NA
data10$Q15A7[data10$Q15A7 == "88"] <- NA
data10$Q15A7[data10$Q15A7 == "99"] <- NA
data10$country[data10$country == "AT"] <- "Austria (AT)"
data10$country[data10$country == "IT"] <- "Italy (IT)"
data10$Vaccinated[data10$Vaccinated == "1_or_2"] <- "1_or_2     "
data10 <- dplyr::rename(data10, Media = Q15A1, Parliament = Q15A2, Healthcare = Q15A3, Government = Q15A4, Science = Q15A5, Schools = Q15A6, Pharma = Q15A7, Country = country)
Box_1 <- melt(data10, id.vars=c("Country", "Vaccinated"))
Box_1$variable <- factor(Box_1$variable, levels = c("Parliament", "Government", "Media", "Pharma", "Schools", "Healthcare","Science"))
Box1 <- ggplot(Box_1, aes(x=variable, y=value, fill=Country)) + geom_boxplot() + scale_fill_manual(values=c("#d7191c", "#2c7bb6"))
Box1
Box_2 <- melt(data10, id.vars=c("Country", "Vaccinated"))
Box_2$variable <- factor(Box_2$variable, levels = c("Parliament", "Government", "Media", "Pharma", "Schools", "Healthcare", "Science"))
Box2 <- ggplot(Box_2, aes(x=variable, y=value, fill=Vaccinated)) + geom_boxplot()

figure <- ggarrange(Box1, Box2, labels = c("A", "B"), ncol = 1, nrow = 2)
figure

### Boxplot Psychosocial
data11 <- select (conjoint2, Q16A1, Q16A2, Q16A3, Q16A4, Q16A5, country, Vaccinated)
data11$Q16A1[data11$Q16A1 == "88"] <- NA
data11$Q16A1[data11$Q16A1 == "99"] <- NA
data11$Q16A2[data11$Q16A2 == "88"] <- NA
data11$Q16A2[data11$Q16A2 == "99"] <- NA
data11$Q16A3[data11$Q16A3 == "88"] <- NA
data11$Q16A3[data11$Q16A3 == "99"] <- NA
data11$Q16A4[data11$Q16A4 == "88"] <- NA
data11$Q16A4[data11$Q16A4 == "99"] <- NA
data11$Q16A5[data11$Q16A5 == "88"] <- NA
data11$Q16A5[data11$Q16A5 == "99"] <- NA 
data11 <- dplyr::rename(data11, hopeful = Q16A1, worried = Q16A2, angry = Q16A3, tired = Q16A4, frustrated = Q16A5, Country = country)
data11$Country[data11$Country == "AT"] <- "Austria (AT)"
data11$Country[data11$Country == "IT"] <- "Italy (IT)"
data11$Vaccinated[data11$Vaccinated == "1_or_2"] <- "1_or_2     "
Box_3 <- melt(data11, id.vars=c("Country", "Vaccinated"))
Box3 <- ggplot(Box_3, aes(x=variable, y=value, fill=Country)) + geom_boxplot() + scale_fill_manual(values=c("#d7191c", "#2c7bb6"))
Box3
Box4 <- ggplot(Box_3, aes(x=variable, y=value, fill=Vaccinated)) + geom_boxplot()
Box4
figure <- ggarrange(Box3, Box4, labels = c("A", "B"), ncol = 1, nrow = 2)
figure

############################## END ##############################
